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Abstract — We discuss applications of nonlinear 
filtering of time series by locally linear phase space 
projections. Noise can be reduced whenever the error 
due to the manifold approximation is smaller than the 
noise in the system. Examples include the real time 
extraction of the fetal electrocardiogram from abdom- 
inal recordings. 

I. INTRODUCTION 

For nonlinear signals with intrinsic instabilities, the 
tasks of noise reduction — or signal separation in gen- 
eral — is difficult to accomplish with the traditional 
spectral approach since the signals may exhibit braod 
band frequency spectra. The theory of nonlinear dy- 
namical systems, or chaos theory, provides alternative 
methods for these purposes based on a phase space 
representation of the data. The theory behind these 
methods is outlined in Ref. in this volume. For the 
derivation of time series methods within the frame- 
work of chaos theory, heavy use has been made of the 
theoretical properties of nonlinear deterministic sys- 
tems. For most methods, this also severely restricts 
the scope of systems they can be applied to. Formally, 
this is also the case for nonlinear filtering procedures 
that exploit the peculiar structure generated by deter- 
ministic phase space dynamics. In this paper we will 
show that if these algorithms are used with care, they 
can also give superior results in situations when pure 
determinism cannot be assumed. The reason is that, 
when represented in a low dimensional phase space, 
also non-deterministic systems may exhibit structures 
suitable for filtering purposes. We will give some ex- 
amples of the latter statement and also discuss prac- 
tical issues that have so far hampered widespread use 
of nonlinear filters. 

II. METHOD 

A scalar time series {s„}, n = 1, . . . , can be unfolded 
in a multi-dimensional effective phase space using time 
delay coordinates s„ — {sn-{m-i)T, ■ • ■ j Sn) (t is a de- 
lay time). If {s„} is a scalar observation of a determin- 



istic dynamical system, it can be shown under certain 
genericity conditions j|, |^ that the reconstructed point 
set is a one-to-one image of the original attractor of the 
dynamical system. We will not assume here that there 
is such an underlying deterministic system. Neverthe- 
less, general serial dependencies among the {s„} will 
cause the delay vectors {s„} to fill the available Tri- 
dimensional space in an inhomogeneous way. Linearly 
correlated Gaussian random variates will for example 
be distributed according to an anisotropic multivariate 
Gaussian distribution. Linear geometric filtering Q 
in phase space seeks to identify the principal direc- 
tions of this distribution and project onto them. The 
present algorithm can be seen as a nonlinear general- 
isation of this approach that takes into account that 
nonlinear signals will form curved structures in delay 
space. In particular, noisy deterministic signals form 
smeared-out lower dimensional manifolds. Nonlinear 
phase space filtering seeks to identify such structures 
and project unto them in order to reduce noise. 

Let us recall the three main steps involved in 
the noise reduction algorithm described in Ref. 
(1) Find a low dimensional approximation to the "at- 
tractor" described by the trajectory {s„}. (2) Project 
each point s„ in the trajectory orthogonally onto the 
approximation to the attractor to produce a cleaned 
vector s„ (3) Convert the sequence of cleaned vec- 
tors s„ back into the scalar time domain to produce a 
cleaned time series s„. 

III. REAL TIME FILTERING 

All of the methods that have been proposed in the 
literature are formulated as a posteriori filters. The 
whole signal has to be available before a cleaned ver- 
sion can be computed, which is then invariably quite 
computer time intensive. One class of methods uses a 
global nonlinear function to represent the dynamics (at 
least approximately). This function has to be deter- 
mined by a delicate fitting procedure and the actual 
filtering scheme (for example Ref. ||^) consists of an 
iterative minimisation procedure. The other class of 



algorithms approximates the dynamics in phase space, 
or phase space geometry, by locally linear mappings. 
Here the tessellation of phase space into small neigh- 
bourhoods is the most time consuming step, along with 
the need to solve a least squares problem in each of 
these neighbourhoods. With fairly low dimensional 
signals, fast neighbour search algorithms (see |6| for 
an overview) are very helpful in this regard. Let us 
introduce some modifications to the locally projective 
noise reduction scheme discussed in Ref. that make 
its use in real time signal processing feasible. (1) The 
data base of local neighbours that is needed to approx- 
imate the dynamics is restricted to points in the past. 
As a side effect, the curvature correction can be carried 
out during the first sweep through the data. (2) The 
number of neighbours required for each point is limited 
to a number that is just sufficient for statistical sta- 
bility. (3) The last modification uses the fact that the 
dynamics is supposed to vary smoothly in phase space. 
The full linearised dynamics is reduced to a collection 
of representative points which is stored together with 
their local linear structure. Consequently, the local lin- 
ear problem has to be solved only for points that are 
not yet well approximated by such a representative. 

Thus, we have turned the procedure into a causal 
filter by restricting neighbour search to points defined 
by measurements made in the the past, k < n. By fur- 
ther limiting the number of neighbours searched for to 
t'^max, we have sped up the formation of the local co- 
variance matrices considerably. Still, as the algorithm 
stands, we have to solve an (m x m) eigenvalue problem 
for each point that is to be processed. A large fraction 
of this work can be avoided on the base of the assump- 
tion that the local linear structure changes smoothly 
over phase space. By making local linear approxima- 
tions we have already assumed smoothness of the un- 
derlying manifold in the Ci sense. In most physical 
systems, the additional assumption of C2 smoothness 
is not less justified. We cannot, however expect that 
the vectors which span the principal directions vary 
slowly from point to point. The reason is that often 
some eigenvalues are nearly degenerate and change in- 
dices from point to point when they are ordered by 
their magnitude. We thus refrain from interpolating 
principal components between phase space points. In- 
stead, we choose a length scale h in phase space which 
is small enough such that the linear subspaces spanned 
by the local principal components can be regarded as 
effectively the same. Now we successively build up a 
data base of representative points for which the local 
points of tangency s*-" , and the local principal direc- 
tions c'^, q ~ 1, . . . , Q have been determined already. 
For each new point s„ that is to be processed, we go 
through this collection of points to determine whether 
a representative is available closer than h. In that 
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Figure 1: Result of nonlinear filtering of an NMR laser 
time series. An enlargement of about one quarter of 
the total linear extent of the attractor is shown. 



case, we use the stored tangent point and principal di- 
rections of the representative in order to perform the 
projections. If not, a neighbourhood is formed around 
s„ in which the eigenvalue problem is solved. The 
point s„ is then included in the list of representatives. 

For real time application, this means that the cor- 
rected value Sn cannot be available before Sn+m-i has 
been measured and processed. Usually, however, this 
delay window is a very short time, at least compared 
to the duration of the recording, and the procedure 
can be regarded as effectively on-line as long as the 
computations necessary to obtain s„ can be carried 
out fast enough. 

IV. RESULTS 

NMR laser data — As a first example, we show 
in Fig. the result of applying the described proce- 
dure to a data set from an NMR laser experiment Q . 
The same data has also been used in Ref. The 
laser is periodically driven and once per driving cy- 
cle the envelope of the laser output is recorded. The 
resulting sampling rate is 91 Hz. At this rate, the non- 
linear noise reduction scheme can be easily carried out 
in real time on a Pentium II processor at 200 MHz. 
Since further iterations can be carried out after the 
time corresponding to one embedding window without 
interfering with the previous steps, we could perform 
up to three iterations on a dual Pentium II worksta- 
tion at 300 MHz in real time. The figure shows the 
result after two iterations. Projections from m = 7 
down to Q = 2 dimensions were used, at least 100 
neighbours were requested at a neighbourhood size of 
2 units. The history was limited to 20000 samples, or 
220 s. This fairly large data base is needed since the 
initial noise level is already small (less than 2% |^) 
and small neighbourhoods are required to avoid the 
dominance of curvature artefacts. 

Magneto-cardiogram data — The actual accel- 
eration resulting from the above modifications strongly 
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a) 


all neighbours 


490 s 


b) 


box assisted neighbour search 


94 s 


c) 


all neighbours in past 


246 s 


d) 


n — k < 5 s 


191 s 


e) 


(d) and C/^ax < 200 


44 s 


f) 


(e) and reuse of representatives 


6 s 



Table 1: Computation time for nonlinear noise reduc- 
tion of 10 s of an MCG recording sampled at 1000 Hz. 
See text for details. 



Figure 2: Result of nonlinear filtering of a MCG time 
series 
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Figure 3: Result of nonlinear filtering of an MCG time 
series 

depends on the situation and it is difficult to give 
general rules and benchmarks. Let us however study 
a realistic example in some detail to illustrate the 
main points. In electrophysiological research, it is 
quite attractive to augment the measurements of elec- 
tric potentials with recordings of the magnetic field 
strength. The latter penetrate intervening tissue much 
more efficiently. Of particular interest are magneto- 
encephalographic (MEG) recordings which allow to ac- 
cess regions of the brain noninvasively which cannot 
be monitored electrically using surface electrodes. In 
cardiology, the magneto- cardiogram (MCG) provides 
additional information to the traditional electrocardio- 
gram (ECG). A particular application is the noninva- 
sive monitoring of the fetal heart which is otherwise 
complicated by shielding of the electric field by in- 
tervening tissue. A common problem with magnetic 
recordings, however, is that the fields are rather fee- 
ble and the measurements have to be carried out in a 
shielded room. Even then, noise remains a major chal- 
lenge for this experimental technique. We will demon- 
strate in the following how nonlinear noise reduction 
could be used for continuous MCG monitoring. 

We use an MCG recording of a normal human sub- 
ject at rest. The data was kindly provided by Carsten 



Sternickel at the University of Bonn. The sampling 
rate was 1000 Hz, which is quite high for cardiac mon- 
itoring. For the signal processing task, any sampling 
rate above about 200 Hz would be sufficient. (Below 
200 Hz, the spike representing the depolarisation of 
the ventricle might not be resolved properly). In or- 
der to cover a significant fraction of one cardiac cycle 
by an embedding window, we chose an embedding with 
delay t = 10 ms in m = 10 dimensions. Neighbour- 
hoods were formed with a radius of 0.06 uncalibrated 
A/D units of the recording, about three times the es- 
timated noise level. We quote computation times for 
the processing of 10 s of MCG on a Pentium processor 
at 133 MHz, determined on a laptop PC running the 
Linux operating system. The timing results are sum- 
marised in Table |^. The data base of representatives 
for (f) was formed by assuming that the local linear 
subspaces are equivalent on length scales of the order 
of 0.06 A/D units. 

Fetal ECG extraction — Let us finally discuss 
the extraction of the fetal electrocardiogram (FECG) 
from non-invasive maternal recordings. Other very 
similar applications include the removal of ECG arte- 
facts from electro-myogram (EMG) recordings (elec- 
tric potentials of muscle) and spike detection in 
electro-encephalogram (EEC) data. 

Fetal ECG extraction can be regarded as a three- 
way filtering problem since we have to assume that a 
maternal abdominal ECG recording consists of three 
main components, the maternal ECG, the fetal ECG, 
and exogenous noise, mostly from action potentials of 
intervening muscle tissue. All three components have 
quite similar broad band power spectra and cannot be 
filtered apart by spectral methods. The fetal compo- 
nent is detectable from as early as the eleventh week of 
pregnancy. After about the twentieth week, the signal 
becomes weaker since the electric potential of the fetal 
heart is shielded by the vernix caseosa forming on the 
skin of the fetus. It appears again towards delivery. In 
Rcfs. |l^, it has been proposed to use a nonlinear 
phase space projection technique for the separation of 
the fetal signal from maternal and noise artefacts. A 
typical example of output of this procedure is shown in 
Fig. H. The assumption made about the nature of the 
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Figure 4: Signal separation by locally linear projec- 
tions in phase space. The original recording (upper 
trace) contains the fetal ECG hidden under noise and 
the large maternal signal. Projection onto the mani- 
fold formed by the maternal ECG (middle) yields fetus 
plus noise, another projection yields a fairly clean fetal 
ECG (lower trace). The data was kindly provided by 
J. F. Hofmeister 



data is that the maternal signal is well approximated 
by a low-dimensional manifold in delay reconstruction 
space. After projection onto this manifold, the mater- 
nal signal is separated from the noisy fetal component. 
Now it is assumed that the fetal ECG is also approx- 
imated by a low-dimensional manifold and the noise 
is removed by projection. Since both manifolds are 
curved, the projections have to be made onto linear 
approximations. For technical details see Refs. §,0. 



V. DISCUSSION 

We have demonstrated that by certain modifications 
to the algorithms described in the literature, nonlinear 
projective noise reduction can be turned into a signal 
processing tool that can in many situations run in real 
time in a data stream. The class of problems that can 
be solved by this approach is much wider than initially 
assumed since strict determinism of the signal is not 
necessary. Signal and noise have to be distinguishable 
by their shape in a reconstructed phase space. 

The algorithms used in this paper and in Ref. I) 
are publicly available as part of the TISEAN software 
project jl^. All the examples were indeed carried out 
with these implementations. 

We thank Led Flepp, Carsten Sternickel, and John 
F. Hofmeister for letting us use their data. This 
work was supported by the SFB 237 of the Deutsche 
Forschungsgemeinschaft . 
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